## Purpose: 

suppressPackageStartupMessages(
  {
    library(tidyverse)
    library(dplyr)
    library(estimatr)
    library(ggplot2)
    library(ggthemes)
    library(texreg)
    library(lfe)
    library(rdrobust)
    library(grid)
    library(gridExtra)
    library(rjson)
    library(tigris)
    library(sf)
    library(tidycensus)
  })


load("crime_la_geo.RData")

### Create rdit

crime_la_ts = 
  crime_la_geo_ts %>% 
  mutate(bg = "fips") %>% 
  mutate(count = 1) %>% 
  group_by(bg, date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_la_ts$blmrv = crime_la_ts$date - min(crime_la_ts$date[crime_la_ts$blm == 1])

rdout_crime_la = with(crime_la_ts, rdrobust(x = blmrv, y = count, p = 1, all = TRUE))

##### Separate terciles to calculate total pop within these #####
## crime
crime_la_geo3 = 
  crime_la_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high")) %>% 
  mutate(bl_3 = ifelse(bl_low3 == 1, "low", "high")) %>% 
  mutate(bg = as.numeric(fips))

crime_la_geo_inc3 = 
  crime_la_geo3 %>%
  group_by(bg) %>% 
  summarise(pop_total = mean(pop_total))

crime_la_geo_inc3 = 
  crime_la_geo3 %>%
  group_by(inc_3) %>% 
  summarise(pop_total_inc3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_la_geo_nw3 = 
  crime_la_geo3 %>%
  group_by(nw_3) %>%
  summarise(pop_total_nw3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_la_geo_bl3 = 
  crime_la_geo3 %>%
  group_by(bl_3) %>%
  summarise(pop_total_bl3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_la_geo3_ts = crime_la_geo3
crime_la_geo_inc3 <- as.data.frame(crime_la_geo_inc3)
crime_la_geo_nw3 <- as.data.frame(crime_la_geo_nw3)
crime_la_geo_bl3 <- as.data.frame(crime_la_geo_bl3)
crime_la_geo3_ts = merge(crime_la_geo3_ts, crime_la_geo_inc3, by = "inc_3")
crime_la_geo3_ts = merge(crime_la_geo3_ts, crime_la_geo_nw3, by = "nw_3")
crime_la_geo3_ts = merge(crime_la_geo3_ts, crime_la_geo_bl3, by = "bl_3")

#### Create plots ####
## crime
# income
crime_la_geo_inclow3 = crime_la_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "A. Crime (Income Lowest Tercile)")

crime_la_geo_inchigh3 = crime_la_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "B. Crime (Income Highest Tercile)")

# Nonwhite pop
crime_la_geo_nwlow3 = crime_la_geo3_ts %>% 
  filter(nw_low3 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "C. Crime (Nonwhite Population Lowest Tercile)")

crime_la_geo_nwhigh3 = crime_la_geo3_ts %>%
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "D. Crime (Nonwhite Population Highest Tercile)")

# Black pop
crime_la_geo_bllow3 = crime_la_geo3_ts %>% 
  filter(bl_low3 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "E. Crime (Black Population Lowest Tercile)")

crime_la_geo_blhigh3 = crime_la_geo3_ts %>%
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "F. Crime (Black Population Highest Tercile)")

crime_la_tercileGrob = arrangeGrob(crime_la_geo_inclow3, crime_la_geo_inchigh3, 
                                   crime_la_geo_nwlow3, crime_la_geo_nwhigh3, 
                                   crime_la_geo_bllow3, crime_la_geo_blhigh3, 
                                   ncol = 3, top=textGrob("Crime by Block Group (LA)"))

ggsave(plot = crime_la_tercileGrob, filename = "crime_geo_la.png", width = 16, height = 8)

#### Analysis: RDs 25, 50, 100 bw
crime_la_geo3_nwlow_25 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-05-03") & date <= as.Date("2020-06-21")) %>% 
  filter(nw_low3 == 1)

crime_la_geo3_nwlow_50 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-04-07") & date <= as.Date("2020-07-16")) %>% 
  filter(nw_low3 == 1)

crime_la_geo3_nwlow_100 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-02-16") & date <= as.Date("2020-09-04")) %>% 
  filter(nw_low3 == 1)

crime_la_geo3_nwhigh_25 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-05-03") & date <= as.Date("2020-06-21")) %>% 
  filter(nw_high3 == 1)

crime_la_geo3_nwhigh_50 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-04-07") & date <= as.Date("2020-07-16")) %>% 
  filter(nw_high3 == 1)

crime_la_geo3_nwhigh_100 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-02-16") & date <= as.Date("2020-09-04")) %>% 
  filter(nw_high3 == 1)

crime_la_geo3_inclow_25 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-05-03") & date <= as.Date("2020-06-21")) %>% 
  filter(inc_low3 == 1)

crime_la_geo3_inclow_50 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-04-07") & date <= as.Date("2020-07-16")) %>% 
  filter(inc_low3 == 1)

crime_la_geo3_inclow_100 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-02-16") & date <= as.Date("2020-09-04")) %>% 
  filter(inc_low3 == 1)

crime_la_geo3_inchigh_25 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-05-03") & date <= as.Date("2020-06-21")) %>% 
  filter(inc_high3 == 1)

crime_la_geo3_inchigh_50 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-04-07") & date <= as.Date("2020-07-16")) %>% 
  filter(inc_high3 == 1)

crime_la_geo3_inchigh_100 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-02-16") & date <= as.Date("2020-09-04")) %>% 
  filter(inc_high3 == 1)

crime_la_geo3_bllow_25 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-05-03") & date <= as.Date("2020-06-21")) %>% 
  filter(bl_low3 == 1)

crime_la_geo3_bllow_50 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-04-07") & date <= as.Date("2020-07-16")) %>% 
  filter(bl_low3 == 1)

crime_la_geo3_bllow_100 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-02-16") & date <= as.Date("2020-09-04")) %>% 
  filter(bl_low3 == 1)

crime_la_geo3_blhigh_25 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-05-03") & date <= as.Date("2020-06-21")) %>% 
  filter(bl_high3 == 1)

crime_la_geo3_blhigh_50 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-04-07") & date <= as.Date("2020-07-16")) %>% 
  filter(bl_high3 == 1)

crime_la_geo3_blhigh_100 = crime_la_geo3_ts %>% 
  filter(date >= as.Date("2020-02-16") & date <= as.Date("2020-09-04")) %>% 
  filter(bl_high3 == 1)

## Separate by outcomes
crime_la_geo3_inclow_ts = crime_la_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_la_geo3_inchigh_ts = crime_la_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_la_geo3_nwlow_ts = crime_la_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_la_geo3_nwhigh_ts = crime_la_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_la_geo3_bllow_ts = crime_la_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_la_geo3_blhigh_ts = crime_la_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_la_geo3_inclow_ts$blmrv = 
  crime_la_geo3_inclow_ts$date - min(crime_la_geo3_inclow_ts$date[crime_la_geo3_inclow_ts$blm == 1])
crime_la_geo3_inchigh_ts$blmrv = 
  crime_la_geo3_inchigh_ts$date - min(crime_la_geo3_inchigh_ts$date[crime_la_geo3_inchigh_ts$blm == 1])
crime_la_geo3_nwlow_ts$blmrv = 
  crime_la_geo3_nwlow_ts$date - min(crime_la_geo3_nwlow_ts$date[crime_la_geo3_nwlow_ts$blm == 1])
crime_la_geo3_nwhigh_ts$blmrv = 
  crime_la_geo3_nwhigh_ts$date - min(crime_la_geo3_nwhigh_ts$date[crime_la_geo3_nwhigh_ts$blm == 1])
crime_la_geo3_bllow_ts$blmrv = 
  crime_la_geo3_bllow_ts$date - min(crime_la_geo3_bllow_ts$date[crime_la_geo3_bllow_ts$blm == 1])
crime_la_geo3_blhigh_ts$blmrv = 
  crime_la_geo3_blhigh_ts$date - min(crime_la_geo3_blhigh_ts$date[crime_la_geo3_blhigh_ts$blm == 1])

datasets_la_geo = list(crime_la_geo3_inclow_ts %>% as.data.frame %>% mutate(count_crimela_inc_low = count) %>% filter(date >= as.Date("2016-01-01")),
                       crime_la_geo3_inchigh_ts %>% as.data.frame %>% mutate(count_crimela_inc_high = count) %>% filter(date >= as.Date("2016-01-01")),
                       crime_la_geo3_nwlow_ts %>% as.data.frame %>% mutate(count_crimela_nw_low = count) %>% filter(date >= as.Date("2016-01-01")),
                       crime_la_geo3_nwhigh_ts %>% as.data.frame %>% mutate(count_crimela_nw_high = count) %>% filter(date >= as.Date("2016-01-01")),
                       crime_la_geo3_bllow_ts %>% as.data.frame %>% mutate(count_crimela_bl_low = count) %>% filter(date >= as.Date("2016-01-01")),
                       crime_la_geo3_blhigh_ts %>% as.data.frame %>% mutate(count_crimela_bl_high = count) %>% filter(date >= as.Date("2016-01-01")))

# outcomes

outcomes = c("count_crimela_inc_low", "count_crimela_inc_high", 
             "count_crimela_nw_low", "count_crimela_nw_high",
             "count_crimela_bl_low", "count_crimela_bl_high")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_la_geo)))

for (i in 1:length(datasets_la_geo)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      mat_bw = matrix(NA, 100, 7)
      
      for(z in c(25, 50, 100)) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_la_geo[[i]][, "blmrv"]),
                   y = datasets_la_geo[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z
        
      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

la_rdit_geo = dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>%
  do.call(rbind.data.frame, .) %>% 
  mutate(Outcome = factor(Outcome,
                          levels = c("count_crimela_inc_low", "count_crimela_inc_high", 
                                     "count_crimela_nw_low", "count_crimela_nw_high",
                                     "count_crimela_bl_low", "count_crimela_bl_high"),
                          labels = c("A. Crime, Low Income",
                                     "B. Crime, High Income",
                                     "C. Crime, Low Percentage Nonwhite",
                                     "D. Crime, High Percentage Nonwhite",
                                     "E. Crime, Low Percentage Black",
                                     "F. Crime, High Percentage Black"
                          )),
         Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                             labels = c("DIM", "Linear", "Quadratic")),
         Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                         labels = c("Triangular",
                                    "Uniform",
                                    "Epanechnikov"))) %>% 
  filter(!is.na(BW))

la_rdit_geo = la_rdit_geo %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                         "Epanechnikov, DIM",
                                                                         "Epanechnikov, Linear",
                                                                         "Epanechnikov, Quadratic")))

la_rdit_geo$BW = as.numeric(la_rdit_geo$BW)

#### Coefficient tests of difference ####
## Create function for s
s <- function(n1, n2, p, s1, s2){sqrt(((n1 - p)*((s1)^2) + (n2- p)*((s2)^2)/(n1 + n2 - 2*p)))}

se <- function(s, SEb1, SEb2, s1, s2){s*sqrt(((SEb1/s1)^2) + (SEb2/s2)^2)}

z <- function(b1, b2, se){(b1-b2)/se}

b_est <- function(b1, b2){(b1-b2)}

s_crimela_25_inc <- s(n1 = 15024, n2 = 4520, p = 2, s1 = 0.0006835595, s2 = 0.0003795159)
se_cla_25_inc <- se(s = s_crimela_25_inc, SEb1 = 0.0006835595, SEb2 = 0.0003795159, s1 = 0.0006835595, s2 = 0.0003795159)
z_cla_25_inc <- z(b1 = -1.094947e-03, b2 = -2.375111e-05, se = se_cla_25_inc)
b_cla_25_inc <- b_est(b1 = -1.094947e-03, b2 = -2.375111e-05)
p_cla_25_inc <- 2*pnorm(abs(z_cla_25_inc), mean = 0, lower.tail = FALSE)

s_crimela_25_nw <- s(n1 = 4432, n2 = 14334, p = 2, s1 = 0.0016008910, s2 = 0.0004973288)
se_cla_25_nw <- se(s = s_crimela_25_nw, SEb1 = 0.0016008910, SEb2 = 0.0004973288, s1 = 0.0016008910, s2 = 0.0004973288)
z_cla_25_nw <- z(b1 = -1.769020e-03, b2 = -4.862888e-04, se = se_cla_25_nw)
b_cla_25_nw <- b_est(b1 = -1.769020e-03, b2 = -4.862888e-04)
p_cla_25_nw <- 2*pnorm(abs(z_cla_25_nw), mean = 0, lower.tail = FALSE)

s_crimela_25_bl <- s(n1 = 3877, n2 = 14387, p = 2, s1 = 0.0013212676, s2 = 0.0004386320)
se_cla_25_bl <- se(s = s_crimela_25_bl, SEb1 = 0.0013212676, SEb2 = 0.0004386320, s1 = 0.0013212676, s2 = 0.0004386320)
z_cla_25_bl <- z(b1 = -1.202441e-03, b2 = -3.932729e-04, se = se_cla_25_bl)
b_cla_25_bl <- b_est(b1 = -1.202441e-03, b2 = -3.932729e-04)
p_cla_25_bl <- 2*pnorm(abs(z_cla_25_bl), mean = 0, lower.tail = FALSE)

s_crimela_50_inc <- s(n1 = 30441, n2 = 9034, p = 2, s1 = 0.0004424702, s2 = 0.0002820733)
se_cla_50_inc <- se(s = s_crimela_50_inc, SEb1 = 0.0004424702, SEb2 = 0.0002820733, s1 = 0.0004424702, s2 = 0.0002820733)
z_cla_50_inc <- z(b1 = -6.618450e-04, b2 = 3.642248e-04, se = se_cla_50_inc)
b_cla_50_inc <- b_est(b1 = -6.618450e-04, b2 = 3.642248e-04)
p_cla_50_inc <- 2*pnorm(abs(z_cla_50_inc), mean = 0, lower.tail = FALSE)

s_crimela_50_nw <- s(n1 = 9038, n2 = 28935, p = 2, s1 = 0.0011004079, s2 = 0.0003231081)
se_cla_50_nw <- se(s = s_crimela_50_nw, SEb1 = 0.0011004079, SEb2 = 0.0003231081, s1 = 0.0011004079, s2 = 0.0003231081)
z_cla_50_nw <- z(b1 = -1.211459e-04, b2 = -2.619865e-04, se = se_cla_50_nw)
b_cla_50_nw <- b_est(b1 = -1.211459e-04, b2 = -2.619865e-04)
p_cla_50_nw <- 2*pnorm(abs(z_cla_50_nw), mean = 0, lower.tail = FALSE)

s_crimela_50_bl <- s(n1 = 8074, n2 = 28986, p = 2, s1 = 0.0008939888, s2 = 0.0003121018)
se_cla_50_bl <- se(s = s_crimela_50_bl, SEb1 = 0.0008939888, SEb2 = 0.0003121018, s1 = 0.0008939888, s2 = 0.0003121018)
z_cla_50_bl <- z(b1 = -4.959026e-04, b2 = 1.198520e-05, se = se_cla_50_bl)
b_cla_50_bl <- b_est(b1 = -4.959026e-04, b2 = 1.198520e-05)
p_cla_50_bl <- 2*pnorm(abs(z_cla_50_bl), mean = 0, lower.tail = FALSE)

s_crimela_100_inc <- s(n1 = 60410, n2 = 17677, p = 2, s1 = 0.0003047255, s2 = 0.0002006155)
se_cla_100_inc <- se(s = s_crimela_100_inc, SEb1 = 0.0003047255, SEb2 = 0.0002006155, s1 = 0.0003047255, s2 = 0.0002006155)
z_cla_100_inc <- z(b1 = 8.737272e-05, b2 = 5.315707e-04, se = se_cla_100_inc)
b_cla_100_inc <- b_est(b1 = 8.737272e-05, b2 = 5.315707e-04)
p_cla_100_inc <- 2*pnorm(abs(z_cla_100_inc), mean = 0, lower.tail = FALSE)

s_crimela_100_nw <- s(n1 = 17581, n2 = 57759, p = 2, s1 = 0.0007729933, s2 = 0.0002172478)
se_cla_100_nw <- se(s = s_crimela_100_nw, SEb1 = 0.0007729933, SEb2 = 0.0002172478, s1 = 0.0007729933, s2 = 0.0002172478)
z_cla_100_nw <- z(b1 = 4.662566e-04, b2 = 1.480816e-04, se = se_cla_100_nw)
b_cla_100_nw <- b_est(b1 = 4.662566e-04, b2 = 1.480816e-04)
p_cla_100_nw <- 2*pnorm(abs(z_cla_100_nw), mean = 0, lower.tail = FALSE)

s_crimela_100_bl <- s(n1 = 15827, n2 = 58019, p = 2, s1 = 0.0006476374, s2 = 0.0002260858)
se_cla_100_bl <- se(s = s_crimela_100_bl, SEb1 = 0.0006476374, SEb2 = 0.0002260858, s1 = 0.0006476374, s2 = 0.0002260858)
z_cla_100_bl <- z(b1 = 6.236811e-04, b2 = 3.630627e-04, se = se_cla_100_bl)
b_cla_100_bl <- b_est(b1 = 6.236811e-04, b2 = 3.630627e-04)
p_cla_100_bl <- 2*pnorm(abs(z_cla_100_bl), mean = 0, lower.tail = FALSE)
